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SUMMARY 

A least-squares method based on the first-order veloci ty-pressure- 
vorticity formulation for the Stokes problem is proposed. This method leads to 
a minimization problem rather than to a saddle-point problem. The choice of 
the combinations of elements is thus not subject to the Ladyzhenskaya-Babuska- 
Brezzi (LBB) condition. Numerical results are given for the optimal rate of 
convergence for equal-order interpolations. 


INTRODUCTION 

The Stokes problem has counterparts in many branches of engineering and 
physics. The Stokes operator is a basic ingredient of more complicated models 
of physical phenomena such as the incompressible Navier-Stokes equations. Any 
good solver of Navier-Stokes equations should at least be able to solve the 
Stokes equations. For this reason, the Stokes equations have been a focal 
point of finite element research for over 20 years. 

For the Stokes flow, the Galerkin mixed method of velocity-pressure formu- 
lation is widely used. However, the mixed method leads to a saddle-point prob- 
lem. Consequently, the combination of velocity and pressure interpolations is 
required to satisfy the Ladyzhenskaya-Babuska-Brezzi (LBB) condition (refs. 1 
to 3), which precludes the application of many seemingly natural pairs of 
velocity and pressure elements. Although various convergent combinations of 
velocity and pressure elements have been developed, most of them are not 
convenient. 

In attempting to overcome this difficulty, Hughes and his colleagues 
(ref. 4) introduced the Petrov-Galerkin method. They put an additional least- 
squares term into the classical Galerkin mixed formulation in order to circum- 
vent the LBB test and to apply an equal-order interpolation. Recently, Hughes 
et al . (ref. 5) have improved this formulation, and a symmetric matrix is 
attained. 
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C.L. Chang, in a paper entitled "A Mixed Finite Element Method for Stokes 
Problem: Acceleration-Pressure Formulation," proposed a method employing 

equal-order interpolations for the Stokes problem. By changing variables and 
transforming the classical velocity-pressure formulation into two first-order 
systems with four and two variables, respectively, the least-squares method is 
applied to the systems sequentially. He has proved the optimal rate of conver- 
gence for both velocity and pressure by using Wendland's approach (ref. 6). 

In this report we develop another least-squares method, a method based on 
the first-order veloci ty-pressure-vorti ci ty formulation. The least-squares 
method leads to a minimization problem and thus is not subject to the restric- 
tion of the LBB condition. It can accommodate equal-order interpolations. 

In order to show the properties of the least-squares finite element method 
for first-order systems, we include an error analysis and numerical experiments 
for a simple one-dimensional problem. 

In the following pages the least-squares method is introduced through a 
simple one-dimensional problem and then an error analysis for this simple prob- 
lem is presented. The veloci ty-pressure-vorti ci ty formulation for the Stokes 
problem is described, followed by an introduction to the least-squares finite 
element method for the Stokes problem. Numerical results showing the optimal 
rate of convergence are then given. 


ONE-DIMENSIONAL PROBLEM 

In order to show the difference between the Galerkin mixed method and the 
least-squares method, we introduce here a simple one-dimensional problem as 


d£u _ , 
d* 2 ‘ f 


i n(0, 1 ) 


(1) 


du 


u(0)=0 ^<1> 


where f 6 L (0,1). It can be rewritten as the following first-order system: 


_ 


dx 


= f 


i n ( 0 , 1 ) 


n du _ n 
p “ dx ~ 0 


u(0) 


i n ( 0 , 1 ) 
p(l) = 0 


( 2 ) 


A weak form of the Galerkin method corresponding to equation (2) is given 
in reference 7. Find u 6 H and p G S, such that 


i 


I 


2 



Vv 6 H 




Vq e s 


( 3 ) 


where H = {v e hUo.I), v(O) = 0}, and S = (q £ L 2 (0,1)}. 

One can construct the finite element method based on equation (3), but 
the interpolation of u and p cannot be chosen independently. To assure 
the existence of the discrete solutions, the following LBB condition must be 
sati sfied: 


V.6S 


h’ 


sup 

» 6H h 





d v 
dx 


dx 




r l 

J 0 


dx j 


(4) 


where Hh, Sh denote the correspond! ng finite element subspace. If we choose 
Hh = set of piecewise linear continuous functions and Sh = set of piecewise 
constant functions, then the rate of convergence is 



(5) 


Now let us consider a least-squares method for equation (2). We construct 
a least-squares functional 


J : X x H -» R 


J(p,u) 



( 6 ) 


where X = {q € H 1 ( 0 , 1 ) ; q(l) = 0}. Taking the variation with respect to p 
and u, and letting 6J = 0 lead to a least-squares weak statement, find 
U = (p,u) G X x H such that 


pl 

J 0 


'dp dq 
.dx dx + 





W = <q,v> G X x H 


(7) 


The corresponding finite element problem is then to find Uh = <Ph.Uh> £ x 
such that 


B(U h ,V h ) = L(V h > VV h = (q h ,v h ) G X h x H h 


( 8 ) 


3 


where 



'V 

dx ’ dx 


dq^ 

dT 


(9) 


( 10 ) 


in which (•,♦) denotes the L 2 inner product and h i s the mesh parameter, 


ERROR ANALYSIS FOR ONE-DIMENSIONAL PROBLEM 
It is easy to verify that 

B(u.v) < c||u|| i • nvfl 1 <n> 

where ||l)|| 2 = ||p|| 2 + ||u|| 2 , in which IH^ denotes the H 1 norm. Thus, 

B(*,*) is continuous on X x H. £}(•,•) is symmetric and the inequalities in 
the Lax-Mi Igram theorem reduce to the single coercivity requirement. There 
exists a constant a > 0 such that for V € X x H 

B(V,V) > a||V||^ (12) 


In fact, 


dx 


q - 


dv 


dx 


B(V,V) = 

in which ||-|| denotes the L 2 norm. Consequently, 

,2 


B(V,V) > 


da 

dx 


From equation (13) we have 


B(V,V) = 


B(V,V) > 


da 

dx 


_ dv 
^ dx 


dv 

dx 


- I) 


da 

dx 


+ nail 2 + llvll 2 + 2 


(fr 0 


(13) 


(14) 

(15) 
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(by integration-by-parts and the boundary conditions in equation ( 2 ), and by 

,2 


the Friedrichs inequality 

B(V,V) > 

From equation (16) 


da 

dx 


, Vv 6 H). So 

2 
ii 

+ v 


dq 

dx + v 


B(V,V> > 

B(V, V) > 

The combination of equations (14) and (17) leads to 


2(B(V,V)) 1/2 > 


dg 

dx 


da 

j + v 
dx 


dg 

dx 


dq 

dx + v 


so that 


4B(V,V) > 

Similarly, the combination of equations (15) and (18) leads to 


2(B(V,V)) 1/2 > 


n dv 
q dx 


dv „ 
dx ' q 


dv 

dx 


(16) 


(17) 

(18) 


(19) 


so that 


4B(V, V) 



( 20 ) 


By combining equations (14), (18), (19), and (20) together, we obtain the 
coerci vi ty 


B<V.V> > {3 (lldllf ♦ ||v|f) > ^ ||V«f (21) 

Therefore, the following theorem about the rate of convergence of the corre- 
sponding finite element solutions can be proved directly (ref. 3). 

Theorem . - Assume f is smooth enough and the finite element interpola- 
tion estimates hold; that is. 


IIP - n h p|| 2 < C p h 2a 

( 22 ) 

II u - n h u|| 2 < c u h 2m 
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then 


IIP - P h ||f + II U - u h || 2 < ah 28 - + h 2m ) (23) 

where denotes the finite element interpolation, Q, and m denote the 
orders of polynomials for p and u, respectively, and C p , C u , and C are 
the constants which do not depend on the mesh size h. 

We may utilize the Aubi n-Ni tsche trick to obtain the optimal L 2 esti- 
mates of the error 


llu - u h || < Ch ktl 

(24) 

lip - p h ll < Ch k+1 


where k = min(&,m). 


From the above procedure, we know that the boundary conditions play an 
important role in guaranteeing the coercivity. Now let us consider the same 
problem with different boundary conditions as follows: 


Case 

(a) 

u(0) = 0, 

u(l) = 0 

Case 

(b) 

u ( 0) = 0, 

p(0) = 0 


(25) 


In case (a), more boundary conditions are given for one variable, and no bound- 
ary condition is prescribed for the other variable. A similar situation is 
found in the Stokes problem, where more boundary conditions are related to 
velocity, and less or no boundary conditions are known for pressure and vortic- 
ity. This makes the verification of solvability of the first-order systems 
nontrivial . 


In case (a), we may follow the same procedure as before to prove the coer- 
civity of the corresponding bilinear form B(V,V). 

In case (b), we still have equations (14) and (15). From equation (14) 
and the Friedrichs inequality ||dq/dx|| 2 > ||q|| 2 , Vq 6 H, we have 


B(V, V) > ||q|| 2 

The combination of equations (15) and (26) leads to 


(26) 


2(B(V,V)) 1/2 > 



_ dv 


dv 

+ 

~ q + dx 

2 

dx 


so that 


♦ 


dv 

dx 


4B(V,V) > 


2 


(27) 



From equation (27) and the Friedrichs inequality for Vv E H, we have 

4B( V , V) > Ml 2 (28) 

By combining equations (14), (26), (27), and (28) together, we obtain the coer- 
civity; that is, for V G H x H 


B< V , V) > L ||V||f 


(29) 


Remark 1 . - The symmetry of the bilinear form B(U,V) and the coercivity 
of B(V,V) guarantee that the matrix for the least-squares finite element 
method is symmetric positive definite. This is an important advantage of the 
least-squares method over the mixed methods. 


Remark 2 . - For the least-squares method, the choice of interpolations for 
p and u is not subject to any restriction as long as 1 > 1, m > 1. How- 
ever, inspection of the second part of equation (2) shows that in order to make 
the residual of this equation equal to zero throughout, one may choose the 
interpolation for u as one order higher than that for p. 


To show the convergence of the least-squares finite element method, we 
performed numerical experiments for three simple, one-dimensional, boundary- 
value problems. We used a uniform mesh containing elements of length h. We 
are interested in the behavior of the error e p = p - Ph and e u = u - Uf-, in 
L? norm for various choices of polynomials of degree Q, for p^ and m for 
u(v Results of numerical experiments are shown in figure 1 and summarized in 
the following computed convergence rates for three model problems. 


Problem 1: u' = p, p' = -x3, u(0) = u(l) = 0 


Order <2 

Order m 

llepll 

ll e u II 

1 

1 

O(h^) 

0(h2) 

2 

1 

0(h|) 

0(h2) 

1 

2 

0(h2) 

0(h2) 

2 

2 

0(h3) 

0(h3) 


Problem 2: u 1 + 3u = p, p ' - 2u = -2x2 + 6x - 2, u(0) = u ( 1 ) = 0 

Problem 3: -(a(x)u)’ = p, p' = f, u(0) = u(l) = 0 

a(x) = - + a( x - x ) 2 
a 0 


f(x) = 2 + 2a( x - x Q ) |tan 

a - 0.5, Xg = 0.5 


«<* - x 0 ) 


+ tan \ocXg)| 


Order Q 

Order m 

l|e p || 

lieu II 

1 

2 

1 

2 

0(h2) 

0(h 3 ) 

0(h2) 
0 ( h 3 ) 
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We obtained the rate of convergence of the method by calculating the norm 
of the error for each h, plotting log ||errorjj versus log h, and calculating 
the slope of this line. All of the computed rates of convergence agree with 
the estimate (eq. (24)). Of course, the error estimate for problem 2 must be 
i proved theoretically. 

| VELOCITY-PRESSURE-VORTICITY FORMULATION OF STOKES EQUATIONS 

Let us consider an application of the least-squares method to the Stokes 
equations. By introducing the vorticity u_= curl u as an auxiliary variable, 
the Stokes equations can be written as 

div u = 0 

v curl u_+ grad p = £ (30) 

u_- curl u = 0 

where u is the velocity vector, p is the pressure, £ is the body force vec- 
tor, and v is the viscosity. We shall consider the following two-dimensional 
problem only: 


3u 3v 
3x + 3y " 


0 


3x 


+ V 



f x 


3p 3w , 
ay * v W f y 


in a 


(31) 


to + 


3u 

3y 


3v 

3x 


= 0 


where to denotes the z-component of u_, (f x ,f y) is the body force vector, and 
Q is a bounded domain in with piecewise smooth boundary r. For accommo- 
dating various combinations of boundary conditions, let {T] J 2 J3>^4,^5) 
denote the sides of T. The unit outward normal vector to r is denoted by 
n, and the tangential vector to r is denoted by t. 

We can write equation (31) in the general form of a first-order system: 

Lu = £ 
or 

3u 3u 

A 1 3x + A 2 3y + = £ (32) 
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where 



‘1 

0 

0 

0‘ 


'0 

1 

0 

O' 

A, = 

0 

0 

0 

0 

1 

0 

0 

-V 

A 2 = 

0 

0 

0 

0 

0 

1 

V 

0 


.0 

-1 

0 

0. 


1 

0 

0 

0. 



'0 0 0 0' 


'0 ■ 


V 

B = 

0 0 0 0 
0 0 0 0 
.0 0 0 1. 

£ = 

f x 

f 

m 

u = 

V 

p 

.CO. 


The boundary conditions should be supplemented to complete the boundary- 
value problem. We may consider the following boundary conditions: 


u,v given on Ti 
Up ,o) given on 
p.ut given on T 3 
p,u given on T 4 
p,u n given on Ts 


(33) 


For example, T] could be the inflow, outflow, or wall boundary; r£, T 3 , T 4 , 
and f 5 could be the free surface, inflow, or outflow boundary. 


The solvability of the boundary-value problem depends on the combination 
of the boundary conditions. For many practical problems, the boundary operator 
does not satisfy the Lopatinski condition (ref. 6). Thus, the proof of solva- 
bility involves a technical difficulty. We address this problem and give the 
error analysis in "An Error Analysis of Least-Squares Finite Element Method of 
Veloci ty-Pressure-Vorti ci ty Formulation for Stokes Problem." (Technical 
Report, Cleveland State University, 1988.) 


LEAST-SQUARES FINITE ELEMENT FORMULATION FOR STOKES PROBLEM 

We first discretize the domain Q as a union of finite elements and then 
introduce an appropriate finite element basis. Let N e denote the number of 
nodes for one element, yj denote the element shape functions, and Ui, v j , p j , 
and coj be the nodal values. Equal-order interpolations are employed, so we 
can write the expansion 


Ne 

u.(x,y) = £ ti(x,y) 
n j-1 3 


U j 

V j 

p j 



(34) 
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We define the least-squares functional 


J : (h’(Q)) -»• R 

J<V - J l|Lu h - f|| 2 - | 

where each element in ^(Q)^ satisfies the boundary conditions. 

By introducing the finite element approximation defined in equation (34) 
into the functional equation (35), we have the minimization problem as follows: 

find u h 6 (hJ(G)) , such that 

J(u h ) < J(y h ) Vv h G (h^O)) (36) 

Obviously, this problem is equivalent to solving the algebraic equations 

KU = F (37) 


3u, 9 Mh 

A 1 37 * A 2 W + f 


( 35 ) 


where the U is the global vector of nodal values, the global matrix K is 
assembled from the element matrices 


K = 


Q. 


(L>V 1 ,Lt 2 


,Lt Ne ) ( L +i * L +2 ’ 


■ L *Ne ) 


dxdy 


(38) 


and the vector F is assembled from the element vector 


in which 


F e - 




<Lt 1 .^2 


.Lt Ne > £ dxdy 


A ■ ,x A l 


♦i .A 


if/.B 


♦ 


l ,x 
0 


i >y 



o 

o 


0 

i 

o 

*1.x 

v *i.y 

*i,y 

-"Vx 

0 

♦ 


(39) 


(40) 


and T denotes the transposition. We observe that the matrix K is symmetric 
positive definite. 


As is customary, we use Gaussian quadrature to evaluate the coefficients 
of K e and F e . The number of "Gauss" points required for the solution is of 
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some importance. Inspection of the fourth part of equation (31) shows that the 
vorticity u and the derivatives of velocity components u and v appear 
simultaneously. When an equal-order interpolation is employed, it is impracti- 
cable to reduce the residual of this equation to zero throughout. For this 
reason, in our numerical experiments we use reduced integration (ref. 8). 


NUMERICAL RESULTS FOR STOKES PROBLEM 

We construct a model problem (ref. 9) corresponding to the polynomial 
divergence-free velocity field 

u(x,y) = x 2 (l - x) 2 (2y -6y 2 + 4y 3 ) 

v(x,y) = y 2 (l - y) 2 (-2x + 6x 2 - 4x 3 ) 

the pressure field 

p(x,y) = x 2 - y 2 

and the vorticity field 

«(x,y) = -x 2 (l - x ) 2 (2 - 12y + 12y 2 ) + y 2 (l - y) 2 (-2 + 12x - 12x 2 ) 

with three groups of boundary conditions shown in figure 2. 

The bilinear element and the eight-node quadratic element were tested and 
uniform meshes were used. We employed one-point Gaussian quadrature for the 
bilinear element, and 2x2 quadrature for the eight-node quadratic element. 

The numerical results of the rate of convergence are shown in figure 3. We 
found that in all tested cases 

ll<gi < ch k+1 ||e p || < ch k *' llejl < ch ktl 

where k is the order of polynomial of shape functions; that is, all variables 

u, v, p, and co converge in L 2 norm at the optimal rate. 


CONCLUSIONS 

In contrast to the Galerkin mixed method based on the veloci ty-pressure- 
vorticity formulation of the Stokes problem (refs. 7 and 10), the least-squares 
method presented here does not depend on the LBB condition. We have shown that 
the least-squares method converges at the optimal rate for equal-order interpo- 
lations. In contrast to the penalty method (refs. 11 and 12) and the Petrov- 
Galerkin method (refs. 4 and 5), the least-squares method does not have any 
added parameters in the scheme. This means that the least-squares method is 
robust. In contrast to the methods based on the stream function-vorti ci ty for- 
mulation (refs. 10 and 12), the least-squares method produces the velocity and 
pressure directly without further numerical calculation. This method can also 
be extended to three-dimensional cases. 

The least-squares methods have been successfully applied to many problems, 
including high-speed flow with strong shocks (B.N. Jiang and G.F. Carey, "A 
Stable Least-Squares Finite Element Method for Nonlinear Hyperbolic Problems, 
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Int. J. Numer. Meth. Fluids, to be published, and ref. 13). It is possible to 
make a general-purpose program based on the first-order system in order to at- 
tack various types of problems. To adapt the method to other problems, it is 
necessary only to modify the subroutines associated with the coefficient matri- 
ces A], A£, and B and the vector function £. The only disadvantage of the 
least-squares method is that each node has more variables. 

The method presented here has already been extended to the incompressible 
Navi er-Stokes equations and will be discussed in another paper. 
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(C) PROBLEM 2. <d) PROBLEM 3. 

FIGURE 1. - COMPUTED CONVERGENCE RATE FOR ONE-DIMENSIONAL PROBLEMS. ORDER OF POLYNOMIALS FOR p AND U. I AND m. 
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FIGURE 2. - BOUNDARY CONDITIONS FOR STOKES PROBLEM 
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